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Abstract 

We consider a boundary value problem in unbounded 2D dou- 
,-C| '_ bly periodic composite with circular inclusions having arbitrary con- 

stant conductivities. By introducing complex potentials, the boundary 
i-G , value problem for the Laplace equation is transformed to a special R- 

C^ I linear BVP for doubly periodic analytic functions. This problem is 

solved with use of the method of functional equations. The i?-linear 
BVP is transformed to a system of functional equations. A new im- 
proved algorithm for solution of the system is proposed. It allows one 
^ ■ not only to compute the average property but to reconstruct the solu- 

tion components (temperature and flux) at an arbitrary point of the 
cn ', composite. Several computational examples are discussed in details 

^— ^ ' demonstrating high efficiency of the method. Indirect estimate of the 

■^ ■ algorithm accuracy has been also provided. 
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1 Introduction 

Heterogeneous media model problems serve the purposes of material science 
studies for the analysis of the various fields and prediction of their properties 
[H [ISl [T71 [26] . Different approaches for study linear inhomogeneous material 
are presented in well-known monographs [H [71 [T31 [ISl [ISl U2\- One of the 
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approaches dealing with composite materials is the so-called homogenization 
method (see [31 [ID])- Mathematical aspects of the higher order homogeniza- 
tion have been developed in |5]. The limiting case for large (close to the 
maximal value) rectangular cross-section cylindrical cavities by means of an 
asymptotic procedure were studied in [3] where explicit analytical expres- 
sions for effective parameters have been also found. Non-local phenomena 
resulting from a high contrast (or anisotropy) of composite structures were 
studied in [Tl|6]. In two- and three-dimensional cases the Rayleigh multipole 
expansions method and its generalizations is effectively used (see, e.g. [22]). 
Various analytical approaches have been discussed in [211 123] . 

Essential progress has been already achieved in the area of numerical 
analysis of composite material properties. Such approach is naturally re- 
stricted to a finite cell (or a few cell - representative element) size in order 
to reduce the computation cost. A vast literature related to this approach 
can be found in [28]. The major advantage of analytical approach is a pos- 
sibility to describe and analyze the material properties by means of explicit 
analytical formulas. This allow one to reveal an influence of the materials 
characteristics (like size, shape, location of components, their material prop- 
erties) on the overall properties of the composite (homogeneous approach) 
[21 m [HI [m [m [201 121] . Recently, the relationship between the effective prop- 
erties in the problem of the heat conduction and elasticity have been revealed 
and effectively exploited [T^ [25]. 

In this paper, we reveal another advantage of the analytical approach 
showing that it is capable to efficiently reconstruct the global and local dis- 
tributions of the physical fields. We consider well-known linear heat con- 
duction problem in 2D unbounded doubly periodic composite with material 
properties independent of the temperature field. The components (inclu- 
sions) are supposed to be disjoint disks formed a doubly periodic structure. 
We consider the steady process governed here by the Laplace equation. We 
will mostly follow by pioneering work [1], but a few important improvements 
will be introduced. First, we slightly change the problem formulation intro- 
ducing more natural conditions at infinity prescribing only average fiux, at 
an arbitrary direction, in contrast to the problem investigated in [1], where 
a special temperature distribution assumed in the direction of the coordi- 
nate system. In the linear formulation, our approach is in fact equivalent to 
periodic conditions for fiux on the boundary of the minimal representative 
cell. 

Although, we treat the problem using the methods developed in [1], reduc- 



ing the corresponding BVP to a system of functional equations with respect 
to doubly periodic analytical functions, we substantially change the algo- 
rithm for the numerical calculations. The algorithm in [^ is mostly oriented 
to find the effective conductivity. In this case, it is sufficient to know val- 
ues of the heat flux in the centers of inclusions only. However, it turns out 
that it may not guarantee the best accuracy when defining values of the flux 
outside the inclusion or reconstructing the temperature distribution. Modi- 
fied algorithm presented in this work allows one to increase an accuracy of 
the numerical computations an any point in the distance from the centers of 
inclusions and to find the temperature field (with accuracy to an arbitrary 
constant). The proposed modification allows us to find the flux distribution 
in an explicit form containing all parameters of the considered model such as 
conductivities of the matrix and inclusions, radii and centers of inclusions, an 
intensity and an angle of the flux. As the previous algorithm, it also relays 
on the values of special Eisenstein functions (|27j). 

The paper is organized as follows. In Section 2 we describe the geometry 
of the considered composites and formulate the mathematical problem basing 
on proper physical assumptions. In Section 3 we briefly overview the auxiliary 
problem stated in [1] , show a connection with the original problem. 

In Section 4 we describe a new algorithm in details and show that both 
components of the solution, the flux and the temperature, can be computed 
in the unique scheme. Finally, numerical calculations are performed and 
discussed in Section 5 to demonstrate algorithm accuracy, robustness and 
effectiveness. 



2 Statement of the problem 

We consider a lattice L which is defined by the two fundamental translation 
vectors "1" and "z" (where z^ = —1) in the complex plane C = M^ (with the 
standard notation z = x + ly). 

Here, the representative cell (see Fig. [1]) will be the square 

<5(o,o) := I 2 = ti + «t2 e C : -- < tp < -, p = 1, 2 I . (1) 

Let £ := IJ {mi + im2} be the set of the lattice points, where mi, m2 G Z. 
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Figure 1: The representative cell Qiofi) within doubly periodic composite 



The cells corresponding to the points of the lattice £ will be denoted by 
Q{mum2) = Q{o,o) +mi + im2:= {z EC: z-mi- im2 E Q{o,o)} ■ (2) 

It is considered the situation when mutually disjoint disks (inclusions) D^ := 
{z E C : I2; — Ofcl < R} with a radius R and boundaries dD^ := {z E 
C : \z — ak\ = R} {k = 1,2, ... ,N) are located inside the cell (5(o,o) and 
periodically repeated in all cells Q{mi,m2)- Let us denote by 



N 



Do := Q(o,o) \ IJ ^'^^ U ^^*^ 



(3) 



vfc=i 



the connected domain obtained by removing of the inclusions from the cell 
g(o,o) (cf. Fig.p. 

Let us consider the problem of determination of the heat flux function of 



a doubly periodic composite material with matrix 

Dmatrix = [J ((-Dq U dQ(^o,0)) + mi + 11712) (4) 



mi,m2 



and inclusions 



N 



Dinc= [j [j {Dk + mi + im2) (5) 

mi ,7712 fc=l 

occupied by materials of conductivities A^ > and Aj > 0, respectively. For 
this purpose, we consider a problem of the determination of the potential of 
the corresponding fields, i.e., a temperature function T = T{x,y) satisfying 
the Laplace equation in each component of the composite material 

AT{z) =0, Z e Dmatrix U Dine, (6) 



1,2,. ..,iV: 



which have to satisfy the following boundary conditions on all dD^, k = 

nt) = n{t), (7) 

m 






nil, 1712 



Here, the vector n = (721,-^2) is the outward unit normal vector to dDk, 
TT- = nT-i- + noS- is the outward normal derivative, and 

an ^ ax '^ oy ' 

Tit) := hm Tiz), T,{t) := lim Tiz). (9) 

z^t,zGDo z^t,z^Dk 

The conditions ©-((H]) form the so-called ideal (perfect) contact conditions. 

The thermal loading for the composite is described weakly by the flux 
given at infinity or more accurately by its intensity A. We assume that the 
flux is directed 6 which does not coincide, in general, with the orientation of 
the periodic cell (see. Fig. [1]). According to the conservation law and the 
ideal contact condition between the different materials, the flux is continu- 
ous in the entire structure. Moreover, as a result of such formulation, the 
temperature, which is also continuous as the results of the ideal transmission 
conditions along the interface between the matrix and inclusions, possesses 
non-zero jumps across any cell. 

In addition, we assume that the heat flux is periodic on y. Thus, 

XmTyix,-) = XmTyix,--] = -AsinO + qi{x), (10) 



where A is the intensity of an external flux. The heat flux is periodic on x, 
consequently, 



A, 



,T^y--,yj =\mTxy-,yj = -Acos9 + q2{y). (11) 



To complement to the average flux conditions at infinity, the latter immedi- 
ately proves that the equalities 

1/2 

qj{Od^ = (12) 

-1/2 

are valid for the unknown functions Qj, {j = 1,2). As a result of flTUl) and 
(ITT|) , the heat flux has a zero mean value along the cell 

^ds = 0. (13) 

on 

^ *^{m-i ,m2) 

From the physics point of view, condition f lT3|) is the consequence of the fact 
that no source (sink) exists in the cells. Moreover, since there is no source 
(sink) inside the composite, i.e., neither in the matrix of the composite, nor 
in any inclusion (the total heat flux through any closed simply connected 
curve is equal to zero), we have 

dT 

-,s = 0. (14) 

We will introduce complex potentials ^{z) and (pk{z) which are analytic 
in Dq and D^, and continuously different iable in the closures of Dq and D^, 
respectively, by using the following relations 

{Re {if{z) + Bz), Z e D,natrix, 
(15) 
■^^^Reifkiz), zeD.nc, 

where B is an unknown constant belong to C. Besides, we assume the the 
real part of (f is doubly periodic in Dq, i.e. 

Re (p{z + 1) — Re ip{z) = 0, Re ip{z + ?) — Re ip{z) = 0. 



Note that in general the imaginary part of if is not doubly periodic in Dq. 

Let us show that ip is single- valued function in Dmatrix- We take a har- 
monic function v in D matrix which is the harmonic conjugate to T. For this 
pair of functions the Cauchy-Riemann equations §^ = |^, ^ = — i^ (or the 

^ -^ ^ ox ay ^ ay ax ^ 

SO called normal-tangent Cauchy-Riemann equations f~ = f^? ^ = ~f^) 
have to be valid. The functions v has the following form: 

{Im {ip{z) + Bz), Z e Dmatrix, 
(16) 
^Im(^fc(z), ^G Anc, 

with the same unknown constant B. 

As it follows from (iT3l) - JHM we can write 

These relations yield that the harmonic function v is single-valued in the 
domain Dmatrix- Therefore, the complex potential i^iz) is single-valued in 

-^matrix ■ 

To determine the flux VT{x,y), we need to obtain derivatives of the 
complex potentials: 

(18) 

Let us rewrite conditions ([7])-([8]) in terms of the complex potentials f{z) 
and ipkiz). Let s be the natural parameter of the curve dD^ and 

d _ d d , . 

ds dx dy 

be the tangent derivative along dD^. Applying the Cauchy-Riemann equa- 
tions, the equality ([H]) can be written as 

X^—{t) = Xi-^{t), \t-ak\=R. (20) 

Integrating the last equality on s, we arrive at the relation 

XmV{t) = XiVk{t) + c, (21) 



where c is an arbitrary constant. We put c = since the imaginary part 
of the function y? is determined up to an additive constant which does not 
impact on the form of T. Using flT6|) . we have 



2A 
lmip{t) = -lm{Bt) + - "-—Iniipkit), \t - ak\ = R. (22) 

'^m + \ 

Using (fTSl) . we are able to write the equahty (jTj) in the following form: 

2A 
Re ^{t) = -Re {Bt) + ^Re ^k- (23) 

Adding the relation ( l23l) and ( l22l) multiplied by ?, and using Re ^pk = '^'° ^'^'^ , Im </)/; 
^^27^, t-ak = j^, we have 



(/p(t) = y^fc(t) - pifkit) - Bt, \t - ak\ = R, (24) 

where p = ^'7^'" . 

Let us now differentiate (EU). First, note that 



W)]' = -(-^^) ip'it), \t-ak\=R. (25) 

This can be easily shown by representing the function if in the form y:>{z) = 

oo 

Y^ ak{z — Ofc)', I2; — Ofcl < -R, and by using the relation t = == + ak on the 

1=0 "''' 

boundary |t — afc| = R. Thus, after differentiating ( 12^ and using ( ITSl) . we 

arrive at the following W-linear conjugation problem ([21]) on each contour 

1^ — fffcl = R, 

R ^' 



m = Mt)+p(- Mt)-B, (26) 

\t — Ofc/ 

with A; = 1,2, ...,A^. 

Remark 2.1 T/ins, the boundary value problem ^-^, fTd\) - fTl\) for har- 
monic functions is reduced to M.-linear conjugation problem ( f^) for analytical 
doubly periodic functions ^, ^1, . . . , tpN- 

We will seek a solution ip{z),ipk{z) of the problem ( !26|) as a sum ^/^(s;) = 
'?/;'^^)(2) + tp'^'^^z), tpkiz) = ipl (z) + ipl {z) of solutions of the following two 
BVPs: 



R 



i^^'\t) = ij'~^\t) + p I ^ ) ^«(t) - B„ (27) 

8 



i^^'\t) = V'f (t) + p (^) ^f (t) - ^B,, (28) 

where i/jI and i/jf, are analytical doubly periodic functions, B = Bi + 1B2. 



3 Formulation of an auxiliary problem 

In this section we briefly overview the auxiliary problem discussed in jl] and 
represent necessary results in convenient for us form. Let T be a solution of 
the boundary value problem ([6])- ([8]) with a constant jump corresponding to 
the external field applied in the x-direction 

f{z + l) = f{z) + l, f{z + i)=f{z). (29) 

The complex potentials ^'•^^(2) and (p), {z) are introduced as follows 

{Re {^'^^\z) + Z), Z e Dmatrix, 
(30) 
^Re^«(2;), zeD,^,. 

Note that ^^^\z) and ^[ (z) are analytic in Dq and D^, and continuously 
differentiable in the closures of Dq and Dk, respectively. Besides, the real 
part of (f^^^ is doubly periodic in Dq, i.e. 

Re^(^n^ + 1)-Re^(^)(^) =0, Re^''^\z + t)-Re^^^\z) = 0. (31) 

In general, the imaginary part of if^^'^ is not doubly periodic in Dq. It turns 
out that they satisfy the following M-linear conjugation boundary value prob- 
lem obtained in [3]: 



^(^)(t) = ^i'\t) - p^[!\t) -t, \t- a,\ = R. (32) 

Differentiating the last equality, we obtain that the boundary value problem 
©-([8]), (!29|) is reduced to the M-linear conjugation boundary value problem 
for analytical doubly periodic functions i(j^^\t(j[ , ■ ■ ■ , V'jv (*^f- S)- 

i^^'\t) = ^^k\t)+p(-^yw^)-i (33) 



with _ 

ri'V flT^ I ^^' ' 1 ^ ^ J-^matrix-i 

and _ 

(35) 

Besides, we mention that when the temperature has a constant jump 
corresponding to the external field applied in the ^/-direction 

f{z+l)=f{z), f{z + i)=f{z)-l, 

the temperature is defined as 

{Re {^'^'^\z) +IZ), Z e Dmatrix, 
(36) 

with corresponding functions ^"^^^ ^^ possess the same properties as the 
functions ^^^^ ^^ , and 



The corresponding M-linear conjugation boundary value problem has a form 

i'^'\t) = i^f\t) + p{j^^W^)-^^ \t-a,\=R. (37) 

The problems fl27j) and fl28l) can be reduced to the problems (l33l) and (1371) 
by the following replacements: 

^^^\z) = B,^^'\z), V^«(z) = B,^^^\z), (38) 

i;^'\z) = B,^('\z), ^f\z) = B,^f\z). (39) 



10 



Remark 3.1 It is easy to verify that the functions '4'^{z) := iip^'^^iz) and 
il'^iz) := zipl (zz) satisfy the following M.-linear conjugation boundary value 
problem 



R ^'- 



^^{t) = i,i{t)+plj-y] V,^(t) + 1, \t-bk\ = R, bk = -iau. (40) 

Note that ip^'^^z) = —iip^{—iz) and ipf\z) = —iil)^{—iz). 

Thus, to find a solution of the problem (126|1 . it is sufficient to find solutions 
'il)''^\z),il))^ {z) and iI)-^{z),iIj^{z) of the problems ([331) ^^^ (HUI) . respectively. 



4 Solution of the problem 

First, let us find the real constants Bi and i?2- 
We introduce further notations 

1 1 

2 2 

I:= j ReV^W(^ + 2y) rfy, J^ := ^ ReV^^(^ + ^i/) d|/. 



In general, the integrals / and /^ differ from zero. As it is shown in Remark 

gj below 

1 

2 

/ lm^^^Ux + -\dx = Q. (41) 

_i 

2 



Taking i?2 = (which corresponds to the problem (p7|) ) and using ([15 
(ITg]) and ([31]), we obtain 

^^92>^ = -A^Im (v^(x + ^)+5)= -A™5iImV^«(x+ ^ 
and 

A™^^^^ = A„Re [ij(^- + tyYB)= A^5i(Re(^(i)(i + zy)+l)). 



11 



Due to (HTj) . integration on [— |, |] the first equality gives 



2 



A., / 5^^rfx = 0. 
dy 



Integrating on [— |, i] the second equahty and applying (fTTj) . we obtain the 

constant i?i: 

_ -Acosg 

""^ - a:;:(Ttt) • ^^'^ 

Similarly, taking Bi = (which corresponds to the problem fl25|l ) and 
using flT^ . f lTS]) and f l5S]l . we obtain 

and 

Using the equality ^^"^'(z) = —i^-^{—iz), we have 



Im U^^^ fx + -))= -Re U^ (- - IX 



I 

T^ + y 



Thus, we get 



Am / — ^ dy = 0. 



2 



aT(x i) 



Integrating on [— i, |] the term A.^ — ^^^^- and applying ( ITOj) . we obtain the 

constant i?2: 

„ —A sin 6^ , , . 

Taking into the account the properties of functions under consideration 
from Sections [2] and [3l the results obtained above and Remark |3. 11 we arrive 
at the following theorem: 

12 



Theorem 4.1 Let T = T(x, y) and T^ = Ti.{x, y) be the solution of the prob- 
lem ([^-([^, OT^) and 07]) . The temperature flux is defined in the following 
form: 

i^rr( \ irr( \ ( ip{z) + B, z = x + ly e Dmatrix, 
dTjx, y) _ ^^ dT{x,y) ^ \ 

^^ ^y \ J^^ Mz), Z = X + tye Dine, 

where 

—A cos 9 A sin 9 



matrix ) 



and 

To find the temperature, it is sufficient to find the functions (yj, (/^i, . . . , yj^y 
(cf. (fT5|) ). These functions can be represented as sums (/^(z) = Lp'^^^{^z) + 
(y9(^)(2;), V9fc(2;) = '~p\,\z) + V^fc (2) of two functions c/?*^^^ and (/^^^^ have to 
satisfy the foUowing BVPs: 



^'-^\t) = ^f{t)-p^'{^{t)-B,t, (45) 



^'~''\t)=^t\t)-p^f{t)-^B,t. (46) 

Analogously to ([38D-(l39D, we have ^^^\z) = Bi^^^\z), ipf\z) = Bi^^^\z) 
and '^'''^\z) = B2^^'^\z), ^f (z) = B2^f (z). It is easy to verify that 

^(2)(z) = ^«(-z^), ^f(^) = ^«(-^^)- 

The functions ^'•^^ and iff, can be found up to an arbitrary constant as 
indefinite integrals of the functions 'ip^^'' and ipl. , respectively (cf. (135|1 ). 
Thus, we arrive at the following statement: 

Theorem 4.2 Let T = T{x,y) and Tk = Tk{x,y) be the solution of the 
problem (0j-(0j, fT^) and fill) . The temperature distribution can be found 
up to an arbitrary constant and is defined in the form fTB) . where 

—A cos 9 A sin 9 

B 



A„(/+l) A„(J^-1) 
13 



-A COS g (1) _ -4 sing ~(i), . 

Now we describe a new algorithm for solution of the problem (1331) . The 
problem fl4UI) can be solved analogously. For convenience, we omit upper 
index in i})^^^ and will write if) below. We shortly describe solvability of the 
problem (15^ using some facts and notation of the paper [1]. 

Notice that we have N contours dD^ and A^ complex conjugation condi- 
tions on each contour dD^ but we need to find A^ + 1 functions ip, tpi, . . . , ip^. 
This means that we need one additional condition to close up the system. 
For this reason we introduce a new doubly periodic function $ which is a 

N 

sectionally analytic in (^(o.o) and in |J Dk and has the zero jumps along each 

fc=i 
dDk, k = 1,2, ... ,N. Such consideration will give an additional condition 

on ip, ipi, . . . , ipN- We will show that $ = 0. 

Let us introduce the sectionally analytic doubly periodic function $ by 

the following formula: 

i^k{z) - P Yl J2* Wm^^yn2,m'lpm{z) - 1, \z - Ofcl < R, 

V'(^) - PY. Y. W^um2,m'ipm{z), Z G Dq, 
m=l mi,m2 

where 



Wmum2,m'ip7n{z) = [ ] Ipm { ^^^^=^^^^ + Clr, 

Z — a^a — nil ~ 'mL2 J \Z — ttra — TTli — 11712 



(48) 



and 

N 



/ J / _, '^mi,m2,m ■ / ^ / ^ ''' m\,m2,rn ' / ^ ''' m\,m2,k- \^) 

m=lmi,r?i2 my^kmi,m2 mi,r?i2 

The "prime" notation in Y means that the summation occurs in all mi 

mi,m2 

and m2 except at [7711,7712) = (0,0). 

Applying A7ialytic Co7iti7iuatio7i Pri7iciple and Liouville's theorem for 
doubly periodic functions, we have that $ = c. 

14 



Let ip and ipk be solutions of the system $(2) = c. Then, in Dq, we have 

^{z) = ^'{z) + c (50) 

with some doubly periodic function ip'. Inserting the last equality in f l35p 
and then in fl30|) . we obtain 

T{z) = Re{^'{z) + cz + z), z e Dq, (51) 

with some function if' which yields c = 0. Thus, we have $(-z) = 0. Writing 
$(2;) = 0, we obtain the following system of linear functional equations 

N 

V^fc(z) = p ^ J2* Wm„m2,mMz) + 1 (52) 

m=l mi ,7712 

which is uniquely solvable with respect to tpk in the space of analytical func- 
tions (for more details cf. |1]). 
The function ijj has the form 

N 
■0(z)=P^ ^ W,ni,m2,mi^m{,z). (53) 

m=l mi,m2 

Let US expand tl>k{.z) into Taylor series 

oo 

^k{z) = Y^ ^ik{z - akY (54) 

1=0 

in order to sum up Wmi,m2,k4'k{z) over all translations mi + 11712. 

The series ^ Wj^k'4'k{z), where j = {mi, 7712) and k is a. fixed number, 
j 
can be represented via the elliptic Eisenstein functions Ei{z) of order / (see 

00 

Yl Wj,kMz) = Y ^ikR'^''''^Ei+2{z - ak). (55) 

j 1=0 



The series Y.' ^j,k'^k{z) := J2 Wj,k^k{z) - [jz^^ ^k (^5= + a^j can 



be written in the form 



Y, ' Wj,kMz) = Y i^ikR''^'^'^(Ji+2{z - ak), (56) 
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where cx/ is the modified Eisenstein function defined by the formula cii^z) := 
Ei{z) — z~K The Eisenstein functions Ei converges absolutely and uniformly 
for / = 3, 4, . . . and conditionally for / = 2 (|27]). 

Thus, we can rewrite the equations ( 13^ and ( 15^ for ip^ and ip as follows: 

N oo oo 

(57) 

N oo 

V'(^) = P E E ^'™ ^'^'""'^ ^'+2(^ - «-)• (58) 

m=l 1=0 

Now we need to find the numerical coefficients ipim of the system fl57|l . 
Note that the equation fl58l) for -j/^ has the same coefficients ipim- Taking a 
partial sum of Taylor series with M first items 



i^k{z) = ipok + i^ik{z - ttk) + ip2k{z - akf H h i^Mk{z - ttk) 



M 



and collecting the coefficients of the like powers of 2 — a^, we obtain the 
formula for definition of ipjk'- 

i^jk = ^i^^^^ , (59) 

where ^j^ is derivative of order j of the function tp^- Then, we get 

•'' m^k 1=0 \ )■ 

M ,, . X. 

+ f E ^^.i^^^'^^H-l)^- ^ .m.rV ^+.+2(0) + /„ (60) 

"'■ 1=0 ^ '' 

Thus, we arrive at the system with N{M + 1) unknown constants ipjk and 
N{M + 1) equations which can be solved numerically. Note that this system 
is obtained for an arbitrary number A^ of inclusions. 



16 



Remark 4.3 Note that doubly periodicity of Ei^2 <md the relations E[j^2 = 

0.5 

-(/ + 2)E;+3, / = 0, 1, 2, . . . imply J Ei+2{x + 0.5z) t/x = /or / = 1, 2, . . . , 

-0.5 
0.5 

while the equality J E2{x + 0.5z) dx = can be obtained by numerical cal- 

-0.5 

0.5 _ 

culation. Therefore, J ip{x + 0.5?) dx = 0. 

-0.5 

Remark 4.4 Note that for finding of the flux distribution (namely, the func- 
tions ip,ipk, ■ ■ ■ 1 "^n) in an explicit form, we change an algorithm of solution 
of the equations ( [J?] ) and 13^) in comparison with the algorithm represented 
in ^. It allows to get more accurate numerical values of the flux in each 
point of considered composite material. 

Remark 4.5 Note that it is possible to find the heat flux distribution in 
each cell Q[mi,m2) with corresponding centers ai + rrii + tm2, a2 + rrii + tm2, 
as + mi + zm2, 04 + mi + zm2. 

Thus in the next section when discussing numerical results we are con- 
centrating only on the computations in the Q(o,o) unit cell. 



5 Numerical results and discussions 

5.1 Accuracy of the computations 

The new algorithm to solve numerically obtained system is realized in Maple 
14 software. As an example, we consider the case when four inclusions are 
situated within one cell, i.e., N = 4. We suppose throughout the compu- 
tations that the heat flux of the fixed intensity A = —1 flows in different 
directions with respect to the main axis. Here the minus sign shows that the 
flux is directed from the right to the left (or from the top to the bottom) 
depending on the angle 9. The conductivity of the matrix is set as A^ = 1, 
while those for the inclusions, Aj, will take different values. 

First we analyze the accuracy of the computations when using the mod- 
ified algorithm described in (IFr|) - ( pUj) . 

We take for our tests a non-symmetrical configuration, with respect to 
Ox-axis, with two inclusions situated close enough to each other as depicted 
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Table 1: Basic parameters computed for M = 4 iterations. Centers of the 
inclusions (the same radii R = 0.145) are defined in (1611) . 



9 


5i 


B2 


1/2 

/ qiiOd^ 

-1/2 


1/2 

/ q2{0d^ 

-1/2 


X, = 0.001 




7r/2 
7r/4 


1.867004 


1.320171 



-1.652566 
-1.168540 


-9.9-10-^ 
2.5-10-12 

-7.0 ■ 10-6 


-3.3 - 10-16 
2.8 - 10-* 
2.0 - 10-* 


Xi = 0.01 




7r/2 
7r/4 


1.840595 


1.301497 



-1.637998 
-1.158240 


-9.3-10-*^ 
5.5-10-12 

-6.6-10-6 


2.2 - 10-16 
2.2 - 10-* 
1.6 -10-* 


Xi = 0.1 






7r/2 
7r/4 


1.6301608 


1.1526977 




-1.511547 
-1.068825 


-5.3-10-6 
1.95 - 10-12 
-3.7-10-6 


-1.1-10-16 
-1.1 -10-* 
-8.1 - 10-9 


A. = 1 




7r/2 
7r/4 


1 



0.7071068 




-1 
-0.7071068 




2.2-10-16 

-1.1-10-16 


2.2 - 10-16 


-1.1-10-16 


A, = 10 




7r/2 
7r/4 


0.66154467 


0.46778272 



-0.61309308 
-0.43352227 


-1.3- 10-6 
-4.8 - 10-13 
-9.1-10-^ 


2.2 - 10-16 
-3.9 - 10-^ 
-2.8-10-^ 


Xi = 100 




7r/2 
7r/4 


0.61045509 


0.43165693 



-0.54238833 
-0.38352646 


-1.7-10-6 
-6.2 - 10-13 
-1.2-10-6 


2.2 - 10-16 
-6.7-10-^ 
-4.7-10-^ 


Ai = 1000 




7r/2 
7r/4 


0.605071574 


0.42785021 



-0.53460516 
-0.37802293 


-1.7-10-6 
-3.1 - 10-12 
-1.2-10-6 



-7.0 - 10-^ 
-5.0 - 10-^ 



on Fig. 121 The centers of the inclusions are situated in the points: 



ai = -0.18 + 0.2Z, 02 = 0.33-0.34z, ag = 0.33 + 0.35z, 04 = -0.18-0.2, 

(61) 
while their radii are the same Ri = R. Note that two inclusions from the 
neighboring cells are situated very close to each other. Thus it is impossible 
to neglect the interactions between them. 



r" 




-0.4 




-0.4 




0.2 



0.4 




Figure 2: The cell Q(ofl)- 



The results of the analysis for different contrast p between the materials 
of matrix and the inclusions are presented in Table 1. The results indicate 
that the maximal error in the computations was on the level of 10~^ or bet- 
ter. As an indirect measures we use here the integrals of the auxiliary func- 
tions Qj showing deviation of the flux within each unit cell from its uniform 
distribution. All these values should be zeros as follows from (fT2l) . Other 
computations below suggest that accuracy of the flux components guarantees 
flve valid digits. 

Remark 5.1 Additional computations performed by us show that for rela- 
tively large values of the radius, an accuracy depends on the positions of the 
inclusions. More symmetric configurations with respect to any axis provide 
better accuracy. For an arbitrary (nonsymmetric) configurations, a higher 
accuracy is achieved for smaller radiuses of inclusions. 
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The flux components Qx and Qy in the center a^ of /c-inclusion 






dTk{ak) 
dx 

dTk{ak) 



2Azr.A 



k^m 



^m + ^k 
'^^k^m 



■ReiJkiak) 



Imipkiak), 



dy Xm. + Afc 

are calculated in accordance with the formula (jBl). The flux components in 
any point of the matrix can be found as 



Q^:^\z) = Xm-Re{4j{z)+B), Q 



(m). 



-A„-Im(^(2:) + 5). 



)irn)l 



)("i), 



We calculate Qx (z) and Qy (z) at the point z = G Dmatrix belonging 
to the matrix. Computations in the Table 2 are given for four consequent 
iterations (M = 1,...,4), and flxed other parameters: 6 = 0, R = 0.135, 



A. 



1, Aj = 100. Note that the accuracy (flve valid digits) predicted from 



the computations in the Table 1 has been conflrmed also from other analysis 
below. Note that the convergence is not monotonous as the consequent terms 
in the series change their signs. 

Table 2: The flux components for different iterations M, while other problem 
parameters are: 9 = 0, R = 0.135, A^ = 1, Aj = 100 and the conflgurations 
of inclusion deflned by fl6T]) . 



M 


Q'i\a,) 


Q'y\a,) 


Q^r\o) 


QT\o) 





1.553864 


-0.005491 


0.779414 


-0.0015248 


1 


1.561173 


-0.005155 


0.769965 


-0.0016577 


2 


1.561668 


-0.005273 


0.764487 


-0.0016911 


3 


1.561562 


-0.005274 


0.764280 


-0.0016916 


4 


1.561566 


-0.005275 


0.764289 


-0.0016913 



One can observe form the table that the computational accuracy is be- 
tween flve or six valid units depending on where the flux is computed. Note 
that the material contract is rather high thus the accuracy is good enough 
inside both the materials (in the matrix and in the inclusions). 



20 



5.2 Average properties of the composite 

The effective conductivity of microscopically isotropic composite materials 
with a small concentration v = NttR^ ^ 1 of periodic inclusions is defined 
by the spherical tensor Ag = Agl where the effective conductivity is computed 
by the classical Clausius-Mossotti formula 

\e = \m-\^^ + 0{u^), Z/-^0, (62) 

1 — pu 

where p = (Aj — \m)/{,\ + A^)- For the history and applications of this 
formula see for example |19] . 

In general case of composites with random non-overlapping inclusions the 
tensor of effective conductivity Ag has a form 

with components obtained in |1] in the case when the temperature has a unit 
jump along Ox-axis (see 029 p ) in the following manner: 

J:-ix:y = xui+'2pu^), (64) 

where parameter \E' depends on the material contrast, p, density of the in- 
clusions, z/, and their positions in the cell: 

^ ^ 1 ^ ~ 

^ = ^(p, z/, aj) = ^ 5^ Ma,). (65) 

i=i 

Using (1571) one can prove that for any fixed p and for any non symmetrical 
distribution of the inclusions, the following estimate is true 

^{p,u,aj) = l + 0{u), u^O. 

Thus, the formula (l62l) coincides with (!64l) with an accuracy of O^u'^) for 
^ = and 6 = 7r/2. Another conclusion which immediately follows from this 
preliminary analysis is that the next asymptotic term in the formula (162!) 
should depend on the positions of the inclusions. 

The effective conductivity in the direction of the coordinate system can 
be found in the form (see Appendix 1): 

X:{e) - iXm = -AX^e-^' + -^ Yl MO, a,) (66) 



N 

k=l 
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for an arbitrary direction of the heat flow propagating in the composite and 
defined by the angle 9. Note that in case of ^ = and A = —1 it coincides 
with the formula flM|) . 

To compare the results obtained on the base of the improved algorithm 
and the classical Clausius-Mossotti formula in values we choose a periodic 
array of four inclusions (A^ = 4) placed symmetrically in the unit cell with 
the centers: ai = -0.25 + 0.25z, as = 0.25 + 0.25z, ag = 0.25 - 0.25z, 
04 = —0.25 — 0.252. Other parameters are: the flux flows along Ox-axis 
{6 = 0), A = —1, inclusions conductivities Aj = 100 and various values 
of the concentration parameter v. We find the effective conductivity using 
formula fl66p by the new algorithms computed the values of function ipi{ai) 
(the imposed symmetry leads to Ag = Ag, A^^ = 0). In Table 3, a comparison 
of results calculated by two methods is presented. 



Table 3: Comparison of the effective conductivities Ag computed by the new 
algorithm and the classic Clausius-Mossotti formula A^^. Here the value 
5\ = (Af ^^ — Ae)/Ae indicates the accuracy of the approximate formula. 



R 


V 


D/d 


AP' 


Ae 




5\ 











1 


1 







0.005 


0.0204 


0.00031 


1.00061607 


1.00061569 


3.84021 


3.8 ■ 10~^ 


0.03 


0.01131 


0.1364 


1.02242010 


1.021928 


3.84362 


4.8 ■ 10"^ 


0.07 


0.06158 


0.3889 


1.12846546 


1.11384 


3.85708 


1.3- 10-2 


0.1 


0.12566 


0.6667 


1.28095768 


1.21935 


3.90145 


5.1 ■ 10-2 


0.12 


0.18095 


0.9231 


1.43123393 


1.30138 


3.96553 


1.0-10-1 



Note that computations for v = 0.0204 exhibits much better accuracy 
that reported in the Tables 1 and 2. This fact has been explained in the 
Remark 15.11 Simultaneously, the classic formula provides the results with 
the relative accuracy of 10"^. On the contrary, the value z/ = 0.18095 in 
our computations stands for rather high inclusions concentration than low. 
Indeed, one can see that in this case the inclusion diameter, D = 0.24, and the 
distance between the inclusions, d = 0.26, are comparable in values. From the 
results presented in the table one can extract the constant in the estimate 
O^u"^) from ( l62l) and the relative error between the classic and the exact 
formulas. As follows from the computation to have the deviation between 
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the formulas less that 1% the distance between the inclusions should be 
larger than their diameters at least two times. Clearly, the results refer to 
the particular chosen configuration of inclusions and the contrast parameter 

P- 

Now we return back to the original asymmetric configuration fl6T|) . The 
new algorithm allows us to defined the tensor of effective conductivity for 
the case too. The components Ag and A^'^ are calculated using the formula 
( 166|) . The components A^ and A^"^ of the tensor (l63l) should be calculated in 
the similar manner by choosing the orthogonal direction of the flux. We use 
for the computations the materials conductivities Aj = 100, Am = 1. Values 
of all components of the tensor Ag as functions on the concentration v are 
presented for different angles in Table 4. 



Table 4: The components of the effective conductivity tensor Ag for the 
configuration of the inclusions given in (I^Tj) for the material constants Aj = 
100 and A™ = 1. 



R 


K 


X-y 


K 


Af 


A:(7r/8) 


A:n^/8) 





1 








1 


0.923880 


0.382683 


0.05 


1.059136 


2.2 -lO^^ 


1.060379 


2.2- 10-^ 


0.978514 


0.405790 


0.11 


1.249047 


4.4 ■ 10-6 


1.272899 


4.4- 10-6 


1.153971 


0.487121 


0.135 


1.348160 


1.1-10-^ 


1.398658 


1.1- 10-5 


1.245542 


0.535253 


0.145 


1.389545 


1.5 ■ 10-5 


1.457632 


1.5 • 10-5 


1.283778 


0.557826 



Thus the composite described by such configuration of the inclusions is 
not isotropic. Fig. |2] gives an idea why this should be in the case as the 
inclusions are placed in the same line if the flux flows vertically while their 
centers are shifted against each other in horizontal direction. Using the 
tensorial rule we additionally checked our computations by computing the 
conductivity of the composite when the flux is directed by the angle 9 = 7r/8. 
The accuracy of the computations using the conductivity tensor Ag and the 
tensorial rule versus the formula (1661) was between 10-^ and 10-^ decreasing 
with the radius growth. 
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5.3 The flux and the temperature distribution 



As an example, we also show the flux distribution inside the cell Q(o,o) for 
different angles and conductivities of inclusions on Fig. [3llH We take for 
calculations the centers of inclusions in the points f l6T]) and the radius R = 
0.145. 



-0.4- 



-0.4 



■ 1 ; 

■ ; / 

■ / / 
^ / / 



■ I 

■ I 
/ / 



/ ^ ^ ^ ^ ^ ' 

J ^ ^ ^ ^ ~ ■ 

/ / ■ 

II 

; I 

I \ 

0.2 0.4 



Figure 3: The flux distribution inside Q(o,o) for \ = 0.01, 9 = 0; n /A. 



According to f fTSj) in order to find the temperature function T = T(x, y), 
one needs to find functions ^p and (^fc, k = 1,. . . ,N. Using flTSj) one can 
recover it, up to arbitrary constant, integrating functions tp and ipk, k = 
1, . . . ,N. To define the constant, we use the boundary conditions fl24p . As 
a result of the integration, Weierstrass zeta-function appears (cf. Appendix 
2). The temperature distribution T{x,y) is presented on Fig. [SHSl The same 
two composite configuration as for the flux distribution are given. Namely, 
we consider the following set of the parameters: A^ = 1, -R = 0.145, ^ = 0; | 
with \i = 100 and A, = 0.01. 

To conclude the paper, we have proposed here the improved algorithm 
to solve the system (HTl) . It allows one to reconstruct the solution and its 
gradient at an arbitrary point of the composite with an extremely high accu- 
racy even for a few numbers of the iteration points M = 4. The constructed 
algorithm is equally effective in finding the average properties of the com- 
posite. We have shown its effectiveness on several examples and discussed 
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Figure 4: The flux distribution inside Q(o,o) &!■ -^j = 100, 6* = 0; 7r/4. 




0.25 




Figure 5: The temperature distribution inside (5(o,o) for -^i = 100, 6^ = 0; 7r/4. 



pecuharities of the computations related to a particular configuration of the 
inclusions. And last but not least, the algorithm constructed in this paper 
may be used to compute practically arbitrary configurations of the composite 
while the respective solution (both its local and average values) may be used 
as benchmarks to test other numerical algorithms. 
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0,M4S 




0.90 




0.38 



Figure 6: The temperature distribution inside Q(o,o) for Aj = 0.01, 9 = 0;j. 

5.4 Appendix 1: Evaluation of the effective conductiv- 
ity. 

The effective conductivity in x-direction can be defined by Green's formula 



as 



TV 



Ki^) = ^rnjj ^ dxdy + XiJ2jJ -Q^ dxdy. 

Do ''=^ Dk 

Using first Green's formula 

{ifjAip + Vif ■ Vtp) dV = (p i) (Vv2 ■ n) dS 



(67) 



L 



(68) 



dU 



with ip = X and (f{x, y) = T in Dq (or (p{x, y) = T^ in the respective domain 
Dk) and taking into account ([6]) we have 



dDo 



dT r dT ^ r dTu 

x——as = (t x——as — > ® x— — as, 

•9^ /aQ(o,o) ^^ t^iJdD, dn 



(69) 



where the curves (9(5(o,o) and dD^ are oriented in the counterclockwise direc- 
tion. 

Here we take into account that the contour integrals along the common 
boundaries dDk annihilate according to (jH])- Last integral can be computed 
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with use of (HU]) and ^ 

ffdT r dT /"^/^ /"^/^ 

// -;-—dxdy = (b x-;-—ds= / x(—Asm9)dx— / x(—Asm9)dx+ 

J J dx Jqq^^^^ dn y_i/2 y_i/2 

Q(o,o) 

I /-i/s I ri/2 

- {-Acose)dy + - {-Acose)dy = -AcosO. 

^ J -1/2 ^ J -1/2 

In the y-direction, we have 

A^(^) = Xn 11 ^ dxdy + XJ2[[^ "^^dy. (70) 

Do ''='^ Dk 

Repeating the same hne of the reasoning with ip = y and ip{x, y) = T in 
Dq (or ip{x, y) = Tfc in the respective domain Dk) with first Green's formula 

m 

dT 

— — dxdy = —A sin 9. (71) 

Q{o,o) 
Thus, we have 

— — dxdy = —A cos —} / / — — dxdy, 
ox ^ J J ox 

Do ''-^ Du 

— — dxdy = —A sin 6' — > / / — — dxdy. 
oy TTiJJ ^y 

Do "-^ Dk 

Then, 

\l - i\l = X^-A cos 9 + iA sin 9) + (A, - A„,) ^ ff (^ - i^\ dxdy. 

^=^ Dk 

Using ( 1T5]) . we have 

\:{9) - i\l{9) = -AX^e-^' + 2A„p^ (f M^) dxdy. 

'='Di 
Due to the mean value theorem for harmonic functions, we finally obtain 
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5.5 Appendix 2 

This section contains the description of the Eisenstein function Ei of first 
order. Properties of the Eisenstein functions of higher orders are described 
in details in [27] . 

The theory of elhptic functions provides the following formula for a lattice 
sum introduced by Rayleigh 



^^ {mi + im2)^"- 



mi,r?i2 



The Eisenstein functions of order p are defines as 

Ep{z) := 2J {z - mi - i'm2)~^ . 



mi,m2 



The function E2 and the Weierstrass function p are related by the identities 

m) 

E^iz) := p{z) + S2, 

where 6*2 = vr for the square array. 

The derivative of the Eisenstein function possesses the following property: 

E'piz) = -p-Ep+iiz). 

Using this relation, we have E[{z) = —E2{z). Thus, integration on [0,2;] 
gives 

Ei{z) = C{z)-7rz, 

where ( is the Weierstrass zeta function, and C'i^) = "pi^)- 
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